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THE  STABILITY  CRITERION  OF  GODUNOV  AND 
RYABENKII  FOR  DIFFERENCE  SCHEMES 

by 

R.  D.  Richtmyer 

ABSTRACT 

This  report  is  based  on  a  lecture  given  at  the  Courant 
Institute  of  Mathematical  Sciences  In  July  19^4,  In  a  special 
series  of  lectures  arranged  by  P.  D.  Lax  on  the  numerical 
treatment  of  boundary  conditions  In  finite  difference  cal- 
culations.  The  stability  criterion  Introduced  by  Godunov  and 
Ryabenkll  as  a  guide  to  the  effect  of  boundary  calculations  on 
stability  Is  described  here  together  with  an  application  made 
by  D.  A.  Quarles,  Jr.  to  shock  fitting  In  fluid  dynamics. 
The  criterion  Is  In  the  spirit  of  the  von  Neumann  condition  In 
that  It  Is  concerned  with  local  normal  modes  of  the  linearized 
difference  equations;  It  Is  formalized  In  terms  of  the  concept, 
defined  by  Godunov  and  Ryabenkll,  of  the  spectrum  of  an  Infinite 
family  of  linear  operators,  and  the  main  theorem  of  the  theory, 
reproduced  In  a  simplified  form  here,  shows  how  to  test  this 
criterion.  In  the  case  of  one  space  variable,  by  the  study  of 
the  local  normal  modes  near  the  boundary,  as  well  as  In  the 
Interior . 


THE  STABILITY  CRITERION  OF  GODONOV  AND 
RYABENKII  FOR  DIFFERENCE  SCHEMES 
by 
R.  D.  Rlchtmyer 

1.    Introduction 

In  the  absence  of  a  rigorous  stability  theory  for 
general  non-linear  difference  schemes,  especially  when  auxiliary 
procedures  such  as  shock  fitting  are  used,  one  generally  turns 
to  simple  criteria,  such  as  that  of  von  Neumann.   Although  that 
criterion  has  been  a  very  useful  guide  to  the  stability  of  those 
parts  of  a  calculation  that  have  to  do  with  the  interior  of  a 
system,  away  from  shocks,  interfaces,  and  the  like,  it  throws 
no  light  on  the  effect  of  boundary  conditions  on  stability. 
Godonov  and  Ryabenkii  have  developed,  in  [1],  [2],  and  [J>] ,    a 
generalization  of  the  von  Neumann  criterion  to  remedy  this 
defect.   Their  generalization,  and  an  application  of  it,  will 
be  described  in  this  report. 

Just  as  for  the  von  Neumann  criterion,  one  takes  the 
first  variation  of,  or  linearizes,  the  difference  equations  and 
looks  for  the  local  normal  modes.   In  the  neighborhood  of  an 
interior  point  of  the  system,  not  on  a  shock,  interface,  on 
other  singularity  the  coefficients  of  the  equations  of  first 
variations  are  almost  constant,  and  the  local  normal  modes  are 
the  Fourier  components.   Near  a  boundary  (internal  or  external). 
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there  are  generally  additional  local  normal  modes  given  by  a 

— >  — > 

superposition  of  one  or  more  terms  of  the  form  exp   ik  x    , 

— > 
where  the  components  of   k   may  have  imaginary  parts  of  such 

sign  that  the  term  decreases  toward  zero  as  one  goes  away  from 
the  boundary,  and  where  the  coefficients  of  the  linear  combin- 
ation are  so  chosen  as  to  satisfy  the  boundary  condition. 
This  may  be  described  as  a  boundary-layer  effect. 

A  very  simple  example  is  provided  by  the  wave  equation, 
which  may  be  written  as  a  system 


(-[]  Su   _  ^  3v      Sv  _  ^u 


and  differenced  in  the  familiar  manner 

,1         1        n         n 
n  +  2      "^    ~   2  ^  1      ~   ^  1 
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for  which  the  von  Neumann  criterion  leads  to  the  Courant- 
Friedrichs-Levy  condition   r  <  1  ,  where   r  =  cAt/Ax  =   To 
To  simplify  the  notation,  we  add  ^  to  all  the  fractional 
indices  and  write  (2)  as 


n   +  1 

n 

-    u  . 

J 

J 

(3) 


r(v5 , 1  -  v;;)  =  0 , 


n  +  1  n  +  1 

^J  +  1  -  ^j  .  1  -  -(^-  +  1  -  ^-  ^  ')  -  0 
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It  is  assumed  that  equations  (3)  hold  for   j  =  0,1,2,...,   and 
that  we  have  in  addition  the  boundary  condition 


>,  ,  n  +  1  ,    n  +  1    „ 

/^  )  Vq     +  qUq     =  0  , 


where   a   is  a  constant.   We  write  [i      for  e      and   A   for 
the  amplification  factor  (see[4]),  and  we  look  for  solutions 
of  (3)  of  the  form 


-      0-  -0 


(A)"^  (ix)'^  , 


where   u*   and  v*   are  constants.   Substitution  into  (3)  and 

cancelling  out  (A)   (i-i)   gives 

(A  -  l)u*  -  r(^L  -  l)v*  =  0 
(6) 

.  (A  -  l)|iv*  -  rA(^L  -  l)u*  =  0 

Equating  to  zero  the  determinant  of  this  pair  of  equations  gives 
the  equation 


o 

(7)  M-^  -  (2  +  IA^lJJ-)  h  +  1  =  0 

r  A 


The  functions  (5)  will  — >  0   as   j  — >  oo   if  and  only  if 
||_l|  <  1  ,  and  it  can  contribute  to  an  instability  only  if 
|a|  >  1  .  The  amplification  factors  of  (3)  for  the  Fourier 
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terms  all  lie  on  the  unit  circle  in  the   A  plane  if  r  <  1  , 
which  is  assumed;  that  Is:   if   |  m- I  =  1  ,  then  (7)  gives 
|a|  =1  .   Therefore,  in  the  region   |a|  >   1    ,    the  solution  of 
(7),  namely 


(8) 


u- 


1  + 


(A-l)' 


2Ar 


2 


+ 


A-1)^  ^  (A-l) 


T 


Ar 


4A  r 


separates  into  two  branches,  on  one  of  which  \\i\    >  1   every- 
where, and  on  the  other   |ix|  <  1  .   Since  we  want   |  i_l  |  <  1  ,  we 
take  the  minus  sign  and  agree  to  remain  on  the  branch  of  (8)  for 
which  the  radical  is   >  0   for,  say,   A  >  1  . 

Substitution  of  (8)  into  one  of  the  equations  (5) 
then  shows  that  we  can  take,  for  example. 


u^ 


(9) 


(A  -  1 


2Ar 


V 


A-1)^  +  4r^A   , 


The  boundary  condition   (4)  requires  that   v*  +  au*  =  0  ;  usin^ 
(9)>  we  can  solve  this  equation  for   A  ,  and  we  find 


(10) 


g  +  rg 
g  +  r 


We  remind  the  reader,  in  passing,  that  since  the  sol- 
ving of  the  equation  v*  +  au*  =  0   involves  a  squaring,  the 
solution  (10)  must  be  substituted  back  into  this  equation  to 
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verify  that  It  is  not  a  spurious  root. 

Clearly,  if  q  >  1  ,  we  have  found  a  local  normal  mode 
that  will   give  an  instability,  where  as  for   |a|  <  1   there  is 
no  such  mode . 

2.    The  Spectrum  of  a  Family  of  Operators; 
The  G-R  Criterion 

The  first  objective  of  Godunov  and  Ryabenkil  was  to 
formalize  what  one  hopes  to  accomplish  by  this  search  for  local 
normal  modes.   To  this  end,  they  define  what  they  call  the 
spectrum  of  a  family  of  difference  operators.   Let 

7    =    |c(At)|0  <  At/    be  a  one-parameter  family  of  linear 
operators  defined  in  a  Banach  space  /^ .       (See  [4].)   The  complex 
number  A   is  said  to  be  a  point  of  the  spectrum  of  the  family  r 
if  there  is  a  sequence   (At).  >  0   and  a  sequence   u^ ^  f^ 
such  that   iK^^gll  =  1  '   ^^"^^^  ->  0  as   ^  ->  oo  ,  and 

(2.1)         ||{c((At)^)  -AI  j  u^ll  ->  0   as  i  ->  ~   . 

If   C(At)   were  an  operator  independent  of  At  ,  A  would  be  a 
point  of  its  spectrum  in  the  ordinary  sense.   In  the  applications 
to  linear  difference  equations,  for  spatially  finite  problems  with 
boundary  conditions,  each  C(At)   can  be  represented  by  a  square 
matrix,  whose  order  increases  indefinitely  as   At  — >  0  .   The 
non-triviality  of  this  definition  is  shown  by  an  example  in  [2], 
arising  from  a  very  simple  partial  differential  equation  in  which 


all  the  matrices   C((At).)   have  Identical  spectra  consisting  of 
a  single  point  (eigenvalue)  whereas   the  spectrum  of  the  family 
y  is  a  closed  disk  in  the   A  plane.   Parter  also  observed  [9] 
that  a  difference  scheme  may  be  such  that  the  spectral  radii  of 
all   C(At)   are  inside  the  unit  circle,  and  yet  be  unstable.   In 
his  example,  the  von  Neumann  condition  is  violated,  as  well  as 
the  Courant-Friedrichs-Lewy  condition. 

For  the  stability  of  the  family  y,    as  defined,  for 
example,  in  [4],  it  is  clearly  necessary  that  its  spectrum  lie 
in  the  (closed)  unit  disk  of  the  complex  plane,  and  the  require- 
ment that  this  be  so  is  the  spectral  criterion  of  stability  of 
Godunov  and  Ryabenkii.   In  the  following,  we  shall  refer  to  it 
as  the  G-R  criterion. 

3'    Comparison  with  the  von  Neumann 
Condition 

If  the  von  Neumann  condition  is  applied  at  any  point 
X   of  continuity  of  the  coefficients  of  the  (linearized)  differ- 
ence equations,  the  limiting  locus,  as   At  — >  0  ,  of  the  eigen- 
values  A   of  the  amplification  matrix  belong  to  the  spectrum 
of  'y    (the  elements   u^   appearing  in  (2.1)  can  be  taken  as 
suitably  chosen  wave  packets).   Therefore,  the  G-R  criterion 
requires   |a|  <_  1  +  0(1)   as   At  ->  0  ,  which  is  weaker  than  the 
requirement   |a|  <  1  +  0(At)  of  the  von  Neumann  condition.   On 
the  other  hand,  if  in  the  neighborhood  of  a  boundary,  there  are 
local  normal  modes  of  the  kind  discussed  in  Section  1,  corres- 


ponding  to  amplification  values  A  with  a  limiting  value   Ap^ 
as   At  — >  0  ,  then   A^   also  belongs  to  the  spectrum  of  T  , 
and  the  G-R  criterion  is  in  this  respect  stronger  than  the 
von  Neumann  conditions.   In  a  sense,  the  G-R  criterion  is  just 
strong  enough  to  take  care  of  the  limiting  behavior,  as 
At  — >  0  ,  of  all  local  normal  modes. 

From  this  point  on,  the  papers  of  Godunov  and  Ryabenkii 
under  review  are  restriced  to  the  case  of  a  single  space  variable 
X  ;  difference  equations  hold  for   0  ^  ^  ^  1  >    ^i^d  there  is  a 
left  boundary  condition  at   x  =  0  and  a  right  boundary  condition 
at   X  =  1  .   The  lemma,  theorem,  and  definitons  quoted  below 
show  how  the  complete  spectrum  of  jr  can  be  found,  in  this  case, 
by  examining  local  normal  modes. 

Comment  1:   The  type  of  local  normal  mode  sought  in 
Section  1  has  to  be  generalized  in  some  cases.   If  \x      is  a 
multiple  root  of  (7)?  for  some  fixed   A  ,  one  may  have  to  con- 
sider, in  addition  to  the  functions  (5)5  which  vary  exponentially 
with  X   as   (i-L )   =  exp  j^ikjAxj  ,  also  functions  that  vary  as 
(ii.)"   times  a  polynomial  in   j  . 

Comment  2:   If   C,  A   and   u   are  the  difference  oper- 
ator, amplification  factor,  and  local  normal  mode  discovered  in 
Section  1,  it  would  seem  that   (  C  -  AlJ  u   equals  zero  rather 
than  merely  approaches  zero  as  in  (2.1).   However,  a  right 
boundary  condition  at   x  =  1   on   J  =  J  would  presumably  not 
be  satisfied  by   u  ,  except  in  the  limit  as   J  — >  00  ;  if 
I  (J.  I  <  1  ,  the  values  of   u  at   J  =  J   decrease  to  zero  as 

I  I J 

Im-I 
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4.    A  Lemma  of  Ryabenkil 

The  following  lemma,  for  equations  with  constant  co- 
efficients, was  given  in  [1]. 

Lemma:   D„   and  D-,   are   p  x  p  matrices  with  real 
elements;  if   det(Dp^  +  i_lD-,  )  4=  0   for  all  [i      such  that   |ii,|  =  1  , 
and  if  ^,    is  a   k-fold  zero  of   detCD^  +  iiD-.)    ,  where  |(j.q!  <  1  , 
then  there  is  a   k-dimensional  subspace   R(m-q)   of  the   p-dimen- 
sional  vector  space  such  that  if  V„  ^  R(iJ.„)  ,  then  the  equations 


(4.1)  ^O^j  +  °1^1+1  =  °  '     J  =  0,1,2, 


have  a  unique  solution  (v.l   such  that   ||V  .  ||  — >  0  as   j  — >  cx)  ; 
furthermore,  there  is  an  a  >  0   s.t.   Ilv  .  11  <  const.   e~  "^   for 
j  =  0,1,2,.. .  . 

Comment  1 :   If  r,  „  ,  with   |ti„|  <  1   is  a   k'-fold 
zero  of   det  {t\Dq   +  D-,  )  ,  there  is  a  similar  space   3(rjp|)   of 
possible  starting  vectors  of  sequences  of  solutions  of  (4.1) 
going  to  zero  as   J  — >  -  oo  , 

Comment  2:   The  proof  given  by  Ryabenkii  is  based 
on  a  canonical  form,  given  by  Gantmacher,  for  pencils  of  matrices; 
an  independent  proof  is  given  below. 

Proof  of  Lemma:   By  hypothesis,   D„  -i-  D,   is  non-singular. 
Therefore,  if 

(^•2)        Dq  +  nD3_  =.  Dq  +  D^  +  (a  -  1)D^  ^  (D^  f  D^)(l  +  (at-  1)C)  , 
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-1, 


where   C  =  ( D„  +  D  )  ~  D,  ,  then  la^   Is  a  k-fold  zero  of 
det  (I  +  (n  -  DO  .   Write 


PJP 


■^  ,  I  +  (H  -  DC  =  P(I  +  (n  -  1)J)P"   , 


where   J  =  J,  +...+  J.   Is  in  Jordan  canonical  form;  then 
M-Q   is  a   k-fold  zero  of   det  (1  +  (^i  -  1)J)  ;  but  this  last 
matrix  has  zeros  everywhere  below  the  main  diagonal  and  ele- 

are  the  eigenvalues  of   C  .   The  determinant  is 


ments   1  +  (p.  -  1)A.   on  the  main  diagonal,  where   A-,,...,',-, 


P 
i=l 


TT.^",  (1  +  (ij.  -  1)AD  ;  precisely   k  of  the  factors  must  be 


proportional  to  la  -  ia.„  ,  hence   1  +  (ti  -  1)A.  =  A.(ij.  -  M-q)  , 
hence   A.  =  1/(1  -  \Xrs)      for   k  values  of   i  .   Consider  one 
of  the  Jordan  boxes   J,,...,  J,   having  1/(1  -  [x.^]      as  eigen- 
value; without  loss  of  generality,  we  may  suppose  it  to  be   J. 
and  to  have  order   q  <_  k  .   Then,  if   v  ,  .  •  .  ,v   are  the  1 
q  columns  of   P  ,   v-,   is  an  eigenvector  of   C   and 
V  , ..,,v    are  generalized  eigenvectors: 


(^.3) 


Cv,  =  T V.  , 

1    1  -  |a.Q   1 


Cv 


Cv 


Vo  +  V.  , 


2    1  -  Hq   2     1 


V   +  V 


q    1  -  iJ-Q   q     q-1 


Therefore,  if  we  call   v^-.  =  0  ,  then  from  [h  .2)  ,    for   t  =  l,...,q. 


12 


(Dq  +  nD^)v^  ^  (Dq  +  D^)  {v^  +  (H  -  1)  (y^ 


V,   +  V,   -, 

\1^     t  t-1 


)}; 


This  Is  an  Identity  in  [o.  ,  so  that 


(4.4) 


Vt  =  (^o^^i^lS 


v^  -  V, 


M-0   ^     ^"^ 


(^•5)     Vt  =  (^0  ^^l^il^^  ^  \-ll 


Call   R(iXq)   the  space  spanned  by  y-^,...,v^    .      For 
V^6R(ia^)  ,  a  solution  of  (4.1)  can  be  obtained  in  the  form 


V  = 
J 


IIt=l  'j  ^t  '   J  =  0,1,...; 


from  (4.4)  and  (4.5), 


q 


0  ^DqV.  +d^v.^^  =  (Dq  +D^)^t=i 


r.. 


,     t         ^0    ^      t  1     V 

a  .  r  +  a  •  ,  1  ^-M'^_ 


,   t     t  , 
+  (-a.  +  a.^^)  v^_^  J  . 


Since   D„  +  D,   is  non-singular  and  the   v,   are  linearly  inde- 
pendent, the  coefficient  of  each  v,   in  the  above  summation 
must  vanish;  i.e. 


t  ~^0   ^  ^  t     1    ^t+1  ^  ^  t+1 

'  •  r  +  a  .  ,  -,  T-  -a  .    +  a  .  ,  -, 

J  ^L^-l     j+1  i-t^-l    J       j+1 


0  ,   t  =  1, . . . ,q-l  , 


'0 


^0 


a^-^  +  a.^    1 


J  M-Q"-^     J"^-^   *^0~"^ 
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Calling   a  .   the  column  vector  whose  components  are   a  .    to 

J  J 

a  .  ,    the  above  becomes 

J 


(^.6) 


^Xo-l 


^J>1 


M- 


0 


^0-1 


where  the  matrix  elements  not  written  are  all  zero.   Since  the 
matrix  on  the  left  is  upper  triangular,  its  reciprocal  is  upper 
triangular  and  has   M-q-1   on  the  diagonal;  therefore  if 'we  write 
a  .  .  =  Ma  .  ,   M  is  upper  triangular  and  has   ll„   on  the  diagonal 

J+1  J  jrr-  D  ^-Q  t. 

The   j —  power  of   M  therefore  has  elements  bounded  by  \i^'^ 
times  a  polynomial  of  degree   q   in   j  ,   hence  bounded  by 
const.  e~  "^    ,    where   ||jl„|  <  e    <  1  . 
Call 


(4.7) 


R 


S  = 


Y_     ©{R(M.i)|  det  (Dq+h^D^)  =  0  ,   ln^l  <  l}  , 
YL     ®{_S(ti.)|  det  (ti.Dq+D^)  =  0  ,  It^.I  <  l}, 


where  the   S(ri,)   are  as  described  in  Comment  1.   Since  the  equa- 
tion  det  (Dq-HiD,  )  =0  is  of  formal  degree  p   and  the  equation 
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det  (r]D„  +  D-,  )   is  obtained  from  it  by  the  substitution  r\    =r  l/\x    , 
the  totality  of  roots  for  which  ,  |  |j.  |  <  1   and   |  t^  |  <  l   (we  have 
assumed  none  =  1  )   are   p   in  number.   (Some  may  be  zero,  of 
course.)   Therefore,   dim  R  +  dim  S  =  p. 

The  solution  of  (4.1)  is  unique  if  the  V.   are  restricted 

J 

to  R  .   If  D-,   is  singular,  then  one  could  add  W  to  V .  -, 
where   D^  W  =  0  and   W€S  .   But  since  the  solutions  of 
D  W.  +  D-,  W  .  ,  =  0  ,  for   W„  €  S  ,  grow  exponentially  with  increas- 
ing  j  ,  the  further  requirement  that   ||V  .  ||  — >  0  as   j  — >  oo 

J 

makes  the  solution  of  (4.1)  unique  for  any  ¥„$  R  . 


5.  '  Formulation  of  the  Problem 

In  [2],  Godunov  and  Ryabenkii  discuss  difference  equa- 
tions of  the  form 


(5'1)        //i\B,  u.,,  =//,aA,  u.,, 
^— (k)   k   j+k   ^^— (k)  k   j+k 


together  with  boundary  conditions  at   J  =  0  and   j  =  J  ;  the 

summations  are  over  a  finite  set  of  values  of  k  ,  the   A,   and 

'  k 

B   are   p  x  p  matrices,  and  the   u.   are   p-component  vectors. 
Such  equations  are  obtained  as  approximations  to  systems  of  linear 
partial  differential  equations  in  x   and   t  . 

The  paper  is  devoted  primarily  to  the  case  in  which 
the   A^  and  B^  depend  on  x   and  Ax  ,  piecewise  continuously 
on  X  and  continuously  on  Ax   for  some  interval 
0  <_  Ax  <_  t  .   The  discussion  below  will  be  restricted  to  the 
much  simpler  case  of  constant  matrices,  although  the  dependence 
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on  Ax   could  be  taken  into  account  easily. 

In  [5],  Godunov  and  Ryabenkli  showed  that  a  system  like 
(5.1)  can  be  reduced  to  the  case  in  which  k    takes  on  only 
two  values,   0  and   1  ,  by  introducing  further  dependent  vari- 
ables, if  necessary.  Just  as  the  scalar  equation 


a.,  +bu.+cu.,T  =0 


can  be  written  as  the  system 


av  .  +  bu  .  +  cu  .  ,  T  =  0  , 
J     J     J+1 


u  .  -  V  .  ,  T  =  0  : 
J     J+1 


also, by  such  reduction,  a  boundary  condition  of  the  form 

\~   o   n+1   X~      n      ,     1    T  J-  o   n+1      n 

/     6,  u,    =  /    a,  u,   can  be  reduced  to  Pr,u„    =  a^u„  . 

^(k)  k  k     ^(j^)  k  k  ^0000 

The  problem  to  be  discussed  is  therefore 

1  "•■  1 

(5.2>       L^__^ \^';tl'L^^^ \ "j,k  (J  =  o,...,j-i) , 


,   ^ ,         Q   n+1       n     „    n+1       n 
(5.3)         V  0  =  "0  ^0  '   Pj  ^  J  =  «J  ''j   ' 


where   p„  ,  a^  ,  Pj  ,  Qj  are  rectangular  matrices  having  p 
columns,  and  where,  in  spite  of  the  notation,   p   and  a., 
do  not  depend  on   J  . 
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In  order  that  the  number  of  equations  be  equal  to  the 
number  of  unknowns,  the  number  of  rows  of   a^   (or  of   P^  )  plus 
the  number  of  rows  of  Qj   (or  of   Pj  )   Is  assumed  equal  to  p 


n 


It  is  assumed  that  when  the   u.   are  given,  (5.2) 

J 

n  -hi 
and    (5.5)    have   a   unique    solution     u    •       (j    =  0,...,J)    ;    that    is, 

there    is   a   matrix      C        of     p(J   +  1)      rows   and   columns    such  that 

{u''j^|  j    =  0,...,JJ    =   Cj    |u|]|  j    =   0,...,jj. 

The  paper  is  concerned  with  the  properties  of  the  family  {c  j 
of  difference  operators. 

A  point   A   of  the  complex  plane  is  called  a  regular 
point  of  the  family  {c  |  if  the  resolvents 


(5.4)         (Cj  -  AI)"^ 


exist  for  all  sufficiently  large   J   and  have  a  bound  that  is 

Independent  of  J  .   All  other  points  clearly  constitute  the 

spectrum  of  the  family  {C,].   .   This  definition  of  the  spectrum 

is  equivalent  to  that  given  in  Section  2. 

In  terms  of  (5.2)  and  (5.3),  the  conditions  on  the 

resolvents  {5-^)    can  be  stated  as  follows:   A   is  a  regular 

point  of  the  system  (5.2),  (5.'^),  if,  for  given   f   -,   ,  0  ,  Y  , 

J  +  2 
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the    equations 


(5.5)  (ABq    -   Aq)    u      +    (AB      -   A    )    u   ^^    =   f      ^ 


(APq    -    Qq)     Uq    =    (t     ,  (APj    -    Qj)     Uj    =    ¥ 


have    a    unique    solution,    for    sufficiently   large      J    ,    and 


(5.6)  ^®^(1)    11^  J  1     ^     Max   I  Max,  .  J|f      ^|    Jl*  ||  ,||y||  |  , 

where      M     does   not   depend   on     f      ^    ,    <D    ,    Y   ,    or      J    . 


3  +  2 


Introducing  matrices 


(5.6|)        D^  =  AB^  -  ^1  '   (i  =  0^1)   ' 
5^  =  Ap^  -  «i  '   (i  =  O'J)   ' 
the  foregoing  condition  on   A   is  simply  that  the  equations 
(5.7)        DqU  +  D-^u  ^^  =  f  _  ^    (j  =  0,...,J-1) 


J-.2 


5qUq  =  0  ,   &jUj  =  ? 


should  he  we 11 -conditioned,  which  Is  defined  to  mean  that  they 
have  a  unique  solution  satisfying  (5.6). 


6.    The  Main  Theorem 

The  main  theorem  of  Godunov  and  Ryabenkll,  for  the 
greatly  simplified  case  of  constant  coefficients.  Is: 

Theorem:   Necessary  and  sufficient  for  (5.7)  to  be 
well-conditioned  for  sufficiently  large   J   is  that 

1)    det  (D„  +  M.D  )  ^   0   for  all  |i.   such  that   |ij,|  =1  . 
2a)    dim  R  =  number  of  rows  of   5„  , 
2b)    dim  S  =  number  of  rows  of   6j  . 
3a)    for  V€R  ,  5qV  =  0   implies   V  =  0  , 
3b)    for  VeS  ,  5jV  =  0   implies   V  =  0  . 
Here   R  and   S   are  the  subspaces  defined  in  connection  with 
the  lemma  of  Ryabenkll.   With  D.   and   6.   interpreted  as  in 
(5-6—),  the  above  conditions  are  equivalent  to  the  G-R  stability 
criterion  for  (5.I), 

Outline  of  proof  of  necessity:   If  any  of  the  above  con- 
ditions is  not  satisfied,  a  set  of  vectors   u.   can  be  found, 

J 

with     max,  ..    \\u.\\   =   1    ,    such  that   the   right    hand   sides    of    (14) 

can  be   made   arbitrarily   small   for   sufficiently  large      J    ,    thus 

contradicting    (5-6).      First,    If    det    (Dq   +  t^Q^^l'    "  "^     with      \[Iq\    =  1    , 

we   take      u.=4-j(l--^)V  [i^'^    ,    where  V     is    such  that 

(D     +  M-qI^i)    V   =  0     and      ||v||    =  1    .      Then    (5.5)    will   give      *   =  0    , 

Y   =  0    ,    f      -j_  =  0(1/J)    .      If      dim  R   >  number   of  rows    of      Oq    , 

then   condition  3a)    of   the   theorem  will   also   fall   and 
Buq^  R    ,     IIuqII    =   1    ,      Sq^O   "  '^    '      ^^^    solution   of 
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DqU   +  D^u   -j^  =  0   gives   *  =  0  ,  f   i  =  0  (  j=0,  .  .  .  ,  J-1 )  ,  and 

J+2 

-a  J 
^  =  0(e    )  .   If   dim  R  <  number  of  rows  of   6„  ,  then  dim  S  > 

number  of  rows  of   6j  ;  hence,  3b)  will  fall. 

Outline  of  proof  of  sufficiency:   For  any   J  ,  the 

system  (5-7)  consists  of  a  finite  number  of  linear  equations 

in  the  same  number  of  unknowns.   Therefore,  if  a  solution  exists 

for  every  choice  of  the  right  hand  sides,  it  is  unique.   One  must 

show  that,  under  each  of  the  following  circumstances,  a  solution 

exists  and  is  suitably  bounded: 


1 )    0  =  0  ,  Y  =  0  ,  f  _   ^  -  0  for   j  :j.  j 


J+2 


2)    0=0,  all   f   1=0 

J+2 


3)   Y  =  0  ,  all  f   1=0 

J+2 


0 


For  case  1),  it  is  first  shown  that  there  are  unique  vectors 
VCR  and  W€  S   such  that 


(6.1)  ^0^  +  ^1^  =  ^ 


J0+  2 


To  do  this,  bases  for   R  and   S  are  chosen  as  in  the  proof 

of  Ryabenkli's  lemma ,  -(^v-,  ,  .  .  .  ,  v  "I  being  the  basis  of   R  and 

fv  -,,...  ,v  \     that  of   S  .   Then  we  write 
\_  r+1'    '    pj 

P 


^    1   =r_  a^  v^  . 


^.t 
Jo+  t    ~  t=l  ' 
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V 


-L 


r  P 

t=l     ^  t=r+l 


Then,  expressing  D^V  by  (4.5)  and  D^W  by  the  corresponding 
equation  concerning  the  subspace   S  ,  where  the  roles  of  Dq 


and  D^   are  reversed,  (6.1)  gives 
1  / 


(6.2) 


where  the  matrix  F   is  upper  triangular  and  has  non-zero  ele- 
ments 1/(1-|-L.)   or   1/(1-^1.)   on  the  main  diagonal,  hence  is 


invertlble . 


Then,   u.   defined  by   u.   =  W  ,  u  .  _^-j^  =  V  and  by 


0 


^0 


(6.3) 


DqUj  +  D-L^j+i  =  0    ^or    j  ^  Jq  , 


u  £  R   for   j  >  Jq+1  ,  Uj^  S    for   j  <_  Jq  , 


satisfy  the  equations  of  the  first  line  of  i^.J);    also   ||u.|| 

J 

decreases  as   j   increases,  for   j  >  j^^  ,  and  as   J   decreases. 


0 


for   J  <  J 


0 


The   u  .   so  contructed  fail  to  satisfy  the  boundary 

J 

conditions,  and  a  convergent  series  of  corrections  is  added; 
the  complete  solution  of  (5.7)  is 


(6.4) 


u  .  +  V  '.  +  W'.  +  V  .  +  w  .  + 

J    J    J    J    J 
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v^   is  chosen  In   R   to  satisfy   5q(uq  +  v^)  =  0  ; 

w^   Is  chosen  in   S   to  satisfy   6j(uj  +  w' )  =  0  ; 

2  p 

Vq   is  chosen  in   R   to  satisfy   Sq^^O  "*"  ^0^  ^  ^  ' 

etc.  (These  choices  are  unique.)   Otherwise  the  v"*:   are  chosen 
in  R  to  satisfy 


Vj  ^  "i^'jrt  =  °  • 


and  the   w  .   are  chosen  in   S   to  satisfy 


0  J     1   j+1 

Then,   [[^"""jll  "^  const.  e~^  ll^'''o'l   ^^^  I|w"'"q||  <   const.  e~^  ||  w""" ..  ||  , 
so  that  the  series  (6.4)  obviously  converges,  for  sufficiently 
large   J  , 

The  solution  for  cases  2)  and  3),  when  only   0   or   Y 
is   =1=  0  ,  is  similar. 

It  remains  to  be  shown  that  the  general  solution  of 
(5.5),  obtained  by  superposing  these  fundamental  solutions, 
satisfies  the  boundedness  condition  (5.6)  for  some  M  .   The  num- 
ber of  terms  in  this  superposition  increases  with   J  ;  nevertheless, 
an  M  can  be  found,  because  greatly  separated  values  of   J   have 
little  influence  on  each  other:   the  fundamental  solution  decays 
exponentially  as   j   departs  from   j   ,  and  the  series 
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00 


^       e-'"^"-^'ol 


=  —  00 


converges . 


7.    Summary  of  the  Procedure 

To  apply  the  G-R  criterion,  we  start  with  the  equations 
In  the  form  (5.I);  It  Is  not  necessary  to  reduce  them  to  the 
standard  form  (5.2).   For  simplicity,  we  shall  restrict  the 
discussion  to  the  case  In  which  the  matrices   A,   and  B,   are 
Independent  of  x  ;  In  general,  they  depend  on  At   or,  what  Is 
the  same  thing,  since  a  relation  Ax  =  g(At)   Is  assumed,  on 
Ax  =  l/j  .   To  show  that  the  G-R  criterion  Is  satisfied,  we 
must  show  that  any  A   with   |a|  >  1   Is  a  regular  point  of  (5'1) 
for  all  sufficiently  small  At  ,  that  Is,  that  the  conditions 
of  the  main  theorem  (Section  6)  are  met  for  all  sufficiently 
large   J  .   (Dp^  and  D-,   depend  parametrlcally  on  A  ,  according 
to  equation  (5.6^-).) 

The  subspaces   R  and   S  ,  which  were  defined  In  con- 
nection with  the  lemma  of  Ryabenkll,  now  depend  on  A   and  on 
At  ,   and  we  write   R(A,At)  ,  S(A,At)  .   However,  It  Is  assumed 
that   A   and  B,   depend  continuously  on  At   In  some  closed 
Interval   0  <_  At  <_  t  ,  so  that  the  above  conditions,  which  have 
to  hold  for  all  sufficiently  small  At  ,  only  need  to  be  verified 
for   At  =  0  ,  and  we  only  need  to  consider  the  spaces   R(A,0) 
and   S(A,0)  .   To  find  R(A,0)  ,  we  substitute  the  expression 
u*(ij.)"^(A)^  for   u"^   Into  (5-1)  with  At  =  0  ,  where   u*   Is 
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a  constant  vector.   This  gives 
(7-1)  M(A,^l)u*  =  0  , 

where   M(A,m.)   is  the  matrix. 


(7.2)  M(A,n)  =  }_         (n)^(AB,  -  A  )  . 

(k)        ^     ^ 


For  (7'1)  to  have  a  non-trivial  solution,  we  need 
(7-3)  det  (M(A,n))  =  0  . 


For  any  complex   A   with   |a|  >  1  ,  let   [i .  ( A )   be  the  roots  of 
(7-3);  we  assume  that   |  p. .  ( A )  |  ^   1  ,  for  otherwise  condition  1) 
of  the  main  theorem  fails,  and  the  equations  are  unstable.   For 
the  given   A  ,  let   [i(^(A)   be  one  of  the  roots  such  that 
|ii„(A)|  <  1  ,  and  let   k  be  its  multiplicity.   We  look  for  vectors 
u*    u*  ,  etc.  satisfying  (7-1): 


(7.4)  M(A,Hq(A))u*  =  0  ; 


if  we  cannot  find  k  linearly  independent  vectors  satisfying 
(7.4)  (this  depends  on  the  rank  of  M(A,^l„(A))  ),  we  look  for 
further  vectors  satisfying 


(7.5)  (M(A,Hq(A)))^u*  =  0  ,   (M(A,i_Lq(A)))5u^  =  0  ,  et 
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and  continue  in  this  way  until  we  have  found   k  linearly 
independent  vectors;  they  span  the  space   R(m.q)   described  in 
the  proof  of  Ryabenkii's  lemma. 

This  process  is  repeated  with  each  of  the  otherroots 
of  (7.3)  for  which   |  i-i .  ( X )  |  <  1  ;  the  total  set  of  resulting 
vectors  span  the  space   R(A,0)  .   We  wish  to  show  that,  for 
this   A  ,  an  arbitrary  linear  combination  V  of  these  basis 
vectors  satisfies  condition  3a)   of  the  main  theorem,  namely 
that  it  can  only  satisfy  the  left  boundary  condition  If  it 
vanishes  identically,  that  is,  that  there  is  no  local  normal 
mode  near  this  boundary  and  having  A   as  its  amplification 
factor.   We  therefore  substitute  this  linear  combination  V 
into  the  boundary  conditions;  this  gives   r   homogeneous  linear 
equations  in  r   unknowns,  where   r  =  dim  R(A,0).   By  condition 
2a)  of  the  theorem,  the  number  of  boundary  conditions  is  equal 
to   r  ;  otherwise  we  should  not  be  doing  all  this  work.   Equating 
to  zero  the  determinant  of  these   r   equations  gives  an  equation 
for  A  ;  if  this  equation  has  a  root  such  that   |a|  >  1  ,  then 
we  have  found  an  unstable  local  normal  mode.   If,  for  each  A 
with  |a|  >  1  ,  there  is  no  such  mode  at  either  boundary,  the 
G-R  criterion  of  stability  is  satisfied. 

Comment  1 :   In  many  practical  cases,  the  matrices 
Ay.     and  B^  can  be  taken  independent  of  At  .   This  is  true, 
for  example,  for  the  Lax-Wendroff  scheme  for  a  system  of  con- 
servation laws,  also  for  all  the  usual  approximations  to  a 
parabolic  equation  without  lower-order  terms,  etc. 
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Comment  2:   Godunov  and  Ryabenkil  also  considered  equa- 
tions with  variable,  but  plecewise  continuous  coefficients.   At 
each  jump  of  the  coefficients,  there  arises  a  similar  possibility 
of  new  local  normal  modes  arising;  hence,  their  form  of  the  main 
theorem,  although  similar  to  the  form  given  in  Section  6,  is 
considerably  more  complicated. 

Comment  3:   Although  the  space   R(A,0)  depends  on   A  , 
its  dimension   r   is  Independent  of   A   for  all   A   such  that 
|a|  >  1  ,  because  the  roots   m,.(A)   of  (TO)  depend  continuously 
on  A   and  the  dimension  of   R(A,0)   can  only  change  by  having 
some  \i        change  from   |^l.|<1   to   |m..|>1;  this  cannot 
happen,  since   |  m.  .  |   is  assumed  never  =  1   for   |a|  >  1  . 

8.    Application  to  Shock  Fitting 

In  problems  of  fluid  flow,  one  type  of  internal  boundary 
is  a  shock  -  a  surface  running  through  the  fluid,  on  which  cer- 
tain flow  quantities  have  Jumps.   The  Rankine-Hugoniot  conditions 
are  equations  connecting  the  magnitudes  of  these  jumps  with  each 
other  and  with  the  shock  velocity,  that  is,  the  velocity  with 
which  the  shock  moves  through  the  fluid.   They  serve  as  boundary 
conditions  connecting  the  solutions  of  the  differential  equations 
on  the  two  sides  of  the  shock.   A  numerical  procedure  that  keeps 
track  of  the  position  of  the  shock  at  each  cycle  of  the  calcula- 
tion and  the  limiting  values  of  the  flow  quantities  on  either 
side  and  applies  the  jump  conditions  as  boundary  condition  for  the 
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finite-difference  calculations  on  either  side  is  called  shock- 
fitting.   In  connection  with  the  development  of  shock-fitting 
procedures  for  multi-dimensional  problems,  the  procedures  that 
have  previously  been  used  for  problems  in  one  space  variable  and 
time  have  lately  been  under  study.   D.A.  Quarles  [5]  has  used 
the  criterion  of  Godunov  and  Ryabenkii  to  study  the  stability  of 
some  of  these  procedures. 

Let   X   denote  a  cartesian  coordinate  in  the  Eulerian 
formulation  of  the  equations  of  motion  of  an  ideal  fluid  (see 
[4],  Chapter  X.),  which,  in  conservation-law  form  (see  [6]),  are 

(8.1)  '    -         ^      U       +^  (m^/p)  +  p  I    =  0   , 

W  \(P  +  e)  m/p^ 

where   p  ,  m  ,  and   e   are  the  mass,  momentum,  and  energy,  res- 
pectively, per  unit  volume;  in  terms  of  the  fluid  velocity   u   and 
specific  internal  energy  c-  ,   we  have   m  =  pu   and   e  =  p  ( c-  +  —  u  ) 
The  pressure   p   is  given  by  an  equation  of  state 

2 

(8.2)  p  =  p(£,p)  =  P(f  -  ^  ,  p)  , 

P    2p 

whereupon  (7-1)  can  be  abbreviated  as 


8-3>         It  "^3?'''"'  =° 


Let  ^  =  ^(t)   be  the   x   coordinate  of  the  shock  at  time   t  , 
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and  let   jJAx    |j  =  l,2,...j    be  the  point  net  used  for  the 
calculation  of  the  smooth  part  of  the  flow  either  side  of  the 
shock,  for  example  by  using  the  Lax-Wendroff  approximation 
(without  viscosity;  see  [6]  or  [?])• 

Let  Pj^     and   ^'^  be  the  values  of   ^(t)   and  of 

^(t)  =  -rrr  ^(t)   at  time   t  =  nAt  .   Let   p^  and   p_   denote 

n 
the  limiting  values  of   p   at  the  shock,  i.e.,  at   x  =  ^   +0 

and   ^'^  -  0  ,  respectively  at   t  =  nAt  ,  and  similarly  for 

p,  u,  and  L  •   The  jump  conditions  are  (we  write  them  for  time 

(n+l)At  ,   on  the  assumption  that  all  quantities  are  known  at 

time   nAt ) : 


,Q    1,  N  /  n+1     ^n-i-1,  n+1     //  nfl     n+l^  /  1 

(8.4)  (u_     -^    )p_    =  /fp_    -  p 


+  /   n+1  ~   n-Hl 


P+     P. 


,Q  r\  /ftJ^+1    n+1  >  n+1 

(0.5)  a  ~   ^+      'P+        =   ^^"^^    > 


(8.6) 


f    n+1  f  n+1  1    /n+1     ,      n+l\/  1  1        \ 

^  -       -  c+      =   2  (^P-     +  p+  ){^  -  ^T+r) 


We  have  assumed  that  the  shock  is  moving  to  the  right  with  respect 
to  the  fluid;  otherwise  the  negative  square  root  must  be  taken. 
For  simplicity,  we  assume  that  the  fluid  ahead  of  the 
shock  is  in  a  known  constant  state,  so  that  quantities  with 
subscript   "+"  are  known  (the  general  case  is  similar).   Then, 
(8.4)  (8.5),  and  (8.6),  together  with  a  kinematical  equation, 
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such  as 


,7)  ^""'  -  ^"  =  1^  ce^^  ^  'e)  , 


and  the  equation  of  state 

,  Q  Q  >  n+1     ,  c  r^+1     1^+1  ^ 

(8.8)  p_    =  p( t  _    .  P_   )  » 

constitute  5  equations  m  the  o  unknowns   u     >  P_    »  P_    ' 
^  n+i  ^  ^n4-   ^  |n+   ^   ^  sixth  equation  is  needed,  and  this  is 

obtained  by  writing  a  difference  equation,  using  (8.1),  con- 
taining differences  between  one  or  more  of  the  flow  quantities 

n  f  1 
at   X  =  ^     and  at  a  nearby  point   x  =  jAx   just  behind  the 

shock,  for  time   (n+l)At  . 

The  question  arises  which  of  the  equations  (8.1) 

should  be  used  for  this  purpose,  or  whether  some  combination 

of  these  equations  should  be.   It  had  been  conjectured  that,  if 

the  equations  are  first  transformed  to  characteristic  form  (see 

[8]) 

(|^+  (u+c)  1^)  p  +  PC  (|^+  (u+c)  1^)  u  =0  , 

(8-9)       (|t--|7)£  -\(|t  ^-It)  P  =0^ 

P 

(|t  +  (-e)  It)  p  -  p^  ^It  ^  (^-)  If)  -  =0  ' 

^      ^   1/2 
Where   c  =  (S^  -^^  +  4^)       is  the  sound  speed,  then  the  first 
p     c     P 

of  (8.9) J  the  one  corresponding  to  the  forward  moving  characteris- 
tic, when  used  as  described  above,  would  lead  to  a  stable  shock- 
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fitting  procedure,  and  the  other  two  would  lead  to  unstable 
ones . 

The  foregoing  sketchlly  described  shock-fitting  pro- 
cedures are  rather  too  complicated  for  stability  analysis, 
owing  to  the  motion  of  the  shock  through  the  point  net,  which 
necessitates  adding  a  new  net  point,  every  few  cycles,  to  the 
part  of  the  net  behind  the  shock  and  somehow  supplying  these 
net  points  with  Interpolated  values  of  the  flow  quantities  with 
which  to  carry  on.   They  are  also  complicated  by  the  fact  that 
the  net  point  Immediately  behind  the  shock  has  to  have   special 
treatment,  since  It  lacks  one  neighboring  net  point. 

For  this  reason  (among  others),  Quarles  considered  a 
modified  method  of  calculation.  In  which  the  Independent  variable 
X   Is  neither  the  Lagnangean  nor  the  Eulerlan  variable,  but  Is 
a  coordinate  In  a  moving  coordinate  system  so  chosen  that  the 
shock  Is  at  a  fixed  value  of  this  coordinate. 

The  details  of  Quarle^s  Investigation  [5]  are  rather 
too  complicated  to  reproduce,  or  even  describe  here,  but  the 
main  conclusions  were  as  follows.   First,  if  the  standard  stab- 
ility condition   for  the  region  behind  the  shock, suitably  modi- 
fied for  the  moving-coordinate  calculation,  is  adhered  to, 
then  the  use  of  the  first  equation  (8.9)  as  described  above  is 
stable,  in  the  sense  that  no  new  local  normal  mode  is  introduced 
with  an  amplification  factor  A   such  that   |a|_>1  .   However,  if 
either  of  the  other  equations  (8.9)  is  used,  a  new  local  normal 
mode  is  found  with   |a|  =1  ,  which  would  suggest  neutral  stability. 
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In  a  numerical  calculation,  on  the  other  hand,  one  cannot  follow 
the  characteristics  exactly;  I.e.,  the  value  of   ( u+c )  ,  or  of  the 
corresponding  quantity  in  the  movlng-coordlnate  calculation,  is 
only  approximately  known.   Quarles  found  that  a  slight  error  of 
the  quantity   (u+c)   would  convert  the  neutral  stability  into 
instability;  hence,  in  this  sense,  either  of  the  second  or  third 
equations  (8.9)  would  lead  to  unstable  shock-fitting  procedures. 

The  foregoing  is  in  accord  with  the  conjectures  mentioned 
earlier,  and  these  conclusions  have  also  been  confirmed  by  numeri- 
cal calculations  made  by  Quarles,  according  to  which  calculations 
based  on  the  first  of  equation  (8.9)  were  stable,  and  the  others 
unstable.   For  details,  see  [5]. 

It  should  perhaps  be  mentioned  that  the  algebraic 
calculations  that  Quarles  had  to  make  were  exceedingly  long 
and  complicated.   They  were  checked  in  collaboration  with  A. 
Lapidus  by  computer  programs,  which  reproduced  the  formal  alge- 
braic manipulations.   It  is  perhaps  reasonable  to  guess  that 
practical  application  of  the  criterion  of  Godunov  and  Ryabenkli 
will  largely  depend  on  computer  aids  to  the  algebraic  manipulation. 
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